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We investigate the chemical dissolution of porous media using a network model in which the system 
is represented as a series of interconnected pipes with the diameter of each segment increasing in 
proportion to the local reactant consumption. Moreover, the topology of the network is allowed to 
change dynamically during the simulation: as the diameters of the eroding pores become comparable 
with the interpore distances, the pores are joined together thus changing the interconnections within 
the network. With this model, we investigate different growth regimes in an evolving porous medium, 
identifying the mechanisms responsible for the emergence of specific patterns. We consider both the 
random and regular network and study the effect of the network geometry on the patterns. Finally, 
we consider practically important problem of finding an optimum flow rate that gives a maximum 
increase in permeability for a given amount of reactant. 

I. INTRODUCTION 

Chemical erosion of a porous medium is a complex process, involving an interplay between flow, transport, reaction 
and geometry evolution. The nonlinear couplings between these processes may lead to the formation of intricate 
dissolution patterns [l|-[3|, the characteristics of which depend strongly on the fluid flow and mineral dissolution rates. 
In particular, in a broad range of physical conditions, long, finger-like channels or "wormholes" are spontaneously 
formed, where the majority of the flow is focused. 

Understanding the details of a dissolution process is of fundamental importance in a variety of geological systems, 
including diagenesis, karst formation, aquifer evolution Q, and melt migration [5|]. It also plays an important role in 
a number of engineering apphcations, such as dam stability @, CO2 sequestration 0, risk assessment of contaminant 

■ migration in groundwater [8] , and stimulation of petroleum reservoirs [H, |l0| . 

Due to the complexity of the problem, analytical results are scarce and limited to either the initial stages of the 

■ process pj]-[l^ or the analysis of simple model systems [15-19]. Thus most studies of the porous media dissolution 
rely on extensive numerical work. A variety of approaches has been followed in modelling, which can be classified into 
three major categories according to the length scales involved. The coarsest description is provided by the Darcy-scale 
models [H, [13, 21] based on continuum equations with effective variables such as dispersion coefficients, Darcy velocity 
and bulk reactant concentrations. On the other side of the spectrum are pore-scale numerical simulations |22l-[2^ where 
the equations for fluid flow, reactant transport and chemical kinetics are solved in an explicitly three-dimensional pore 
space. Naturally, the pore-scale models are much more accurate than the Darcy-scale ones but also highly expensive 
computationally, and thus limited by the system sizes that they can represent. Finally, somewhere in between these 
two levels of description are the network models [9|, lID., 25], which model fluid flow and dissolution in a network of 
interconnected pipes, with the diameter of each network segment or pipe increased in proportion to the local reactant 
consumption. This is also the approach followed in the present work. 

In constructing the model, we are largely following the ideas of Hoefner and Fogler [9[ but with a few important 
modifications. First, we take into account a potential limiting role of the mass transfer of reactants to the pore 
surface. In that way, we obtain the description of the system in terms of two dimensionless parameters: the effective 
Damkohler number, Dae//, relating the reaction rate to the mean fluid velocity in the pores and another dimensionless 
number, G, measuring the extent to which the dissolution rate is hindered by diffusive transport of reactant across 
the aperture. 

Another important novel element that we introduce is to allow for dynamically changing topology of the network 
as the dissolution proceeds. As the diameters of the eroding pores become comparable to the interpore distances, 
the pores are joined together, thus changing the interconnections within the network. This allows for a more realistic 
representation of the evolving topography of the dissolving porous medium. 

Although the original motivation for the construction of the Hoefner and Fogler network model [9I] was the analysis 
of acidization experiments in limestone cores, another natural application of the model is the simulation of early stages 
of karst formation in the limestone bedding planes which separate the individual strata in the rock. As elucidated by 
Ewers [2^ and Dreybrodt ^ in a first approximation one may regard a bedding plane as a two-dimensional porous 
medium with an average pore size comparable to the grain size of the confining rock, thus 2d network models seems 
to be particularly suited for simulations of such a system. 

The paper is organized as follows: In Sec. [IT] we describe the network model used to represent the evolving porous 
medium. Results of the numerical simulations are given in Sees. IIIIIIVl First, in Sec. IIIII we analyze the form of the 
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dissolution patterns as a function of flow velocity and surface reaction rate. In particular, we identify the regime 
in which the hierarchical structure of dissolution channels is formed and study their length distribution. Then, in 
Sec. [iVlwe consider the dissolution of regular lattices. Finally, in Sec. [Vlwe consider a problem of optimal injection 
rate, important in petroleum reservoir stimulation, so as to achieve the maximum increase in permeability for a given 
amount of reactant ^S, dO|, HO, IHI, iH . We flnish with a summary of our results and conclusions. 



II. THEORETICAL MODEL 

The theoretical model used in this work is an extension of the Hoefner and Fogler scheme [ol], who have pioneered 
the apphcation of pore network models to study dissolution processes. Pore-network modelling has been extensively 
used in simulating single- and multiphase flow in porous media (see e.g. books by Dullien [29| and Sahimi [30[ for a 
detailed description of different models and the discussion of their predictive abilities). However, the applications of 
pore models to dissolution processes remain scarce [ol, |lo|, HH, HH . 

In the model, the porous medium is represented as a triangular network of cylindrical tubes. The points where the 
tubes meet are referred to as nodes. The nodes can either be placed on a regular hexagonal lattice of lattice constant 
Iq or randomly, as illustrated in Fig. [H In the latter case, a random displacement, uniformly sampled from a value 
between — 0.4/o and 0.4/o, is added to the lattice nodes. 

Each pore is a tube with an initial diameter do, which gets enlarged during the dissolution process. A reactive fluid 
is injected into the network through a set of inlet nodes^ where the pressure Pinit) is imposed, and leaves the system 
through outlet nodes where the pressure is kept at zero, pout = 0. Depending on the physical situation we intend 
to model, the inlet pressure is either kept constant, or it is being adjusted in each time step so as to keep the total 
flow through the system constant. The concentration of reactant in the incoming fluid is kept at a constant value of 
c = Cin- The inlet and outlet nodes can either be positioned at separate points within the network or placed along 
its two opposite boundaries, as illustrated in Fig. O Additionally, periodic boundary conditions are applied along the 
lateral direction. 



Figure 1: Triangular lattice with regular or random node configuration. 



A. Flow problem 



Fluid flow in the pores is described by the Hagen-Poiseuille equation for the volumetric flux 
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where {Pj — Pi) denotes pressure drop along the pore joining node i with node j, qij is the volumetric flux in this 
pore, dij and Uj are its diameter and length respectively and fi is the dynamic viscosity of the fluid. At each node, 
we also have a continuity condition 
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Figure 2: A schematic representation of a network with Une inputs and outputs. 



dx 




qc(x+dx) 



Figure 3: A schematic view of a pore and reactions at its surface. 



where the sum is over all the nodes connected by a pore with node j. The resulting system of sparse linear equations 
can then be solved for pressure values at the nodes. 

B. Dissolution of a single pore 

Let us now consider one of the pores comprising the network and study the evolution of its diameter as a function 
of time. We assume that the walls of the pore are dissolving according to the linear rate law 

J f kCyj . 

Here Jr is the reactive flux (the number of absorbed reactant molecules per unit area and unit time), k is the surface 
reaction rate and - the reactant concentration at the wall. The precise meaning of these variables depends on the 
particular reaction of interest, for example for carbonate dissolution by a strong acid, the respective reaction is 

CaCOs + H+ ^ Ca^+ + HCO3 . (4) 

In this case c is a concentration of ions and Jr describes the flux of hydrogen ions as illustrated in Fig. [3l However, 
when calcite is dissolved by aqueous CO2 at pH values similar to those of natural groundwater, the relevant variable 
is rather the calcium ion undersaturation f32l, |33| . Note that in general the full chemistry of the carbonate dissolution 
is more compl ex |32.] . but a single reactant description is a reasonable approximation in a broad range of chemical 
conditions fiol. [33j. 

Additionally, we are taking into account the fact that, in order to take part in the reaction, the reactant first needs 
to diffuse from the bulk towards the pore surface. The diffusive flux can be expressed in terms of the difference 
between the surface concentration, c^^, and the bulk concentration, c by using a mass-transfer coefficient or Sherwood 
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number |34j , 

■J.„^^ic-c., (5) 

Here, the fluid bulk mean concentration of the react ant, c, also known as the mixing cup concentration is defined in 
terms of flow-weighted average of the concentration field over the cross-section of the pore: 



d 

1 



c{r)v{r)27rrdr. (6) 



In general, the Sherwood number, Sh, depends in a complicated way on reaction rate at the pore surface {k) but 
the variation is relatively small [3J, [35|, bounded by two asymptotic limits: high reaction rates (transport limit), 
Sh = 4.364, and low reaction rates (reaction limit), Sh = 3.656. In the numerical calculations we approximate Sh 
by a constant value Sh = 4. Additionally, we neglect entrance effects, which otherwise make the Sherwood number 
dependent on the distance from the inlet of the pore, x. 

By equating the reactive and diffusive fluxes Jr = Jdiff we can relate the wall concentration to the bulk one 

a, 

" 1 + kd/DSh- 

which can then be used to express the reactive flux in terms of the bulk concentration in the pore 

Jr = keffC, (8) 

where the effective reaction rate is given by 

= l + M/DSh' 

As seen from above, the relative importance of reactive and diffusive effects is described by the ratio 

g{d) = kd/DSh (10) 

In sufficiently narrow pores the diffusion is fast and able to keep the concentration field almost uniform across the 
pore diameter. In such a case, characterized by ^ <C 1, the dissolution becomes reaction limited and k^ff ~ k. In the 
other limiting case, for wide pores and/or fast reactions, the reaction rate becomes hindered by diffusive transport of 
reactant across the aperture. As ^ 1, the dissolution can become entirely diffusion limited with keff ~ DSh/d. 

Next, the concentration of the reactant along the pore can be obtained from the mass balance equation. Neglecting 
axial diffusion, this leads to 

dc 

^dx ~^^^e//C, (11) 

where x is the coordinate along the axis of the pore. For a pore of a constant diameter, this can be solved to yield 

c{x) = coe 5 , (12) 

with Co denoting the concentration of reactant at the inlet of a pore. 
Finally, the erosion rate of the pore surfaces is described by 

dt{d/2) = ^c, (13) 

where Cgoi is the concentration of soluble material and u accounts for the stoichiometry of the reaction. The total 
volume of a mineral dissolved from the walls of a pore of diameter d over time At is then 

AV = ll^h^U^ t c{x)dx = Mq^{l - e-'''^^«'/9). (14) 

VCsol Jo CsolV 
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In order to keep the model tractable, we assume that the pores are broadening uniformly along their length. Thus 
the above volume change corresponds to the enlargement of the pore diameter by 

Ad=^ = ^^(1 - e-'^'^^^//'/'). (15) 
Trdl TTdl Csoi^ 

As seen from the above, the decay of the concentration along the length of the pore is determined by a function 

M,) = ^^ = -^, (16) 

q 1 + g{d) 

which relates the reaction rate to the rate of convective transport. Using / and g one can rewrite formulae (p!5)) and 
(p!2|) in a simple form 

do (1+^)/^ ^' ^ ^ 

c{x) = coe-f^^K (18) 



and 



where the time has been scaled by 



and 



i = 2k-ft/do, (19) 



7 = ^ (20) 
is the acid capacity number or volume of solid dissolved by a unit volume of react ant. 



C. Dimensionless parameters characterizing the evolution in the network 

The functions / and g characterizing the dissolution in the pores depend on time through d and q. Additionally, 
their values can vary across the sample, since the pores have different lengths and diameters and carry different flows. 
In the following, while characterizing the effects of flow rate and surface reaction rate on the dissolution patterns, we 
will construct the phase diagram of the patterns in terms of characteristic values of these parameters at t = 

Trdoklo/qin . . 

^^eff = fe^^ (21) 

^ ShD 

and 

G = kdo/DSh, (22) 

where Iq is the lattice constant of the underlying hexagonal network, do is the initial pore diameter, whereas qin is a 
characteristic flow in the inlet pores, i.e. 

Qin = Q/Ninlet-, 

where Niniet is the number of inlet pores. The first of these parameters. Dag//, can be interpreted as an effective 
Damkohler number relating the reaction rate to the mean fluid velocity in the pores, whereas G, as has already been 
mentioned, measures the extent to which the dissolution rate is hindered by diffusive transport of reactant across the 
aperture. The parameter G plays a similar role to the Thiele modulus in chemical engineering [36| which measures 
the relative importance of diffusion and reaction. 
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Figure 4: Change of the network topology associated with pore merging. When the sum of the diameters of the pores AC and 
BC becomes larger than 2/o the points A and B are joined and the two pores are replaced by a single pore with a diameter 
equal to dAC + dsc- 

D. Implementation of the model 

The above-described model is implemented numerically in the following way: 

• Pressure in each node and flow through each tube are calculated from equations ^ and (|2]) using MUltifrontal 
Massively Parallel Solver (MUMPS) m^^. 

• The concentration field in each pore, starting from the inlet ones, is then obtained by a repetitive use of Eq. (p!8|) 
under the assumption that at the pore intersections the flow is divided accordingly to the total volume flux 
through each tube 

• Finally, the diameters of the pores are updated according to Eq. (p!5]) . 

Additionally, we introduce the possibility of merging the pores when their diameters become comparable to the 
interpore distances. To be more specific, when di -\-dj becomes larger than 2/o, we replace both pores by a single pore 
of diameter d = di -\- dj^ as illustrated in Fig. HJ Note that in this way we assure that the total reactive surface area 
before and after merging is the same, the condition which is crucial to guarantee that the overall react ant balance will 
not change as an effect of merging. On the other hand, the total volume of the pores is not conserved during merging 
and the overall hydrodynamic resistance decreases. This, however, has only a minor effect on the system, since at the 
moment of merging the diameters of the pores are large and thus their resistances are almost negligible. 

An important parameter controlling the merging process is the ratio of the initial pore diameter do to the lattice 
constant, Iq. Since the condition for merging is di + dj > 2/o, in the systems with larger do/lo the pores will merge 
faster. In real rocks, the values of pore aspect ratio span a rather broad range, depending on the type of porosity 
present in the rock structure [39]: whether it is a network of interconnected microcracks ('crack porosity' with aspect 
ratios as low as 0.1%) or intergrain pore space ('equant porosity' with aspect ratios of 0.1 — 1). 

III. CHARACTERIZATION OF DISSOLUTION PATTERNS 

Fig. [5] illustrates typical dissolution patterns for 200 x 200 random lattice over a range of different flow and reaction 
rates. There is one inflow in the system and three outflows situated at equal distances from the inlet, but not 
symmetrically. The pores are initially uniform in diameters. The widths of the lines representing the pores in the 
figure are proportional to their diameters. For the sake of clarity, we plot only the pores which have reached d = 3do 
at a given time. The frames in the figure correspond to the 'breakthrough' moment when the dissolution reaches the 
outlet of the system, that is at least one of the outlet pores have broaden 3 times. Unless otherwise stated the value 
of pore aspect ratio was taken to be do/lo = 0.025. 

The phase diagram is plotted in terms of dimensionless parameters Dag// and G, defined in the preceding section. 
For small Dag//, due to the low reaction rate and large flow, the reagent spreads evenly throughout many parallel 
pores along the main flow path and erodes them almost uniformly, which results in a diffuse, turnip-shaped pattern. 
This can be rationalized by noting that the characteristic penetration length of the reactant is, according to Eq. (fT8|) 

lp = lo/Bd.eff (23) 

thus for Dag// ^ 10~^ — 10~^ most of the pores along the main flow path will be invaded almost instantaneously. 



Figure 5: The patterns emerging in the dissolution of 200 x 200 lattice with one inlet and three outlets (marked by circles) 
in a range of different Dae// and G numbers. The effective Damkohler number is the main factor responsible for the way the 
react ant is distributed throughout the system, whereas G controls the amount of dissolution in the pores near the inlet. The 
simulations were performed under the fixed pressure drop conditions. The frames in the figure correspond to the breakthrough 
moment, and the corresponding dimensionless time, t — 2kjt/do, is shown for each frame. 




Figure 6: Dissolution patterns obtained by Ewers in the experiments on the dissolution of blocks of Paris plaster [2^. The 
inlet and outlet of the fluid are marked by circles. 

However, as Dae// is increased, the penetration length is reduced and only the pores in the close vicinity of the 
inlet are invaded at the beginning of the dissolution process. The reagent is then consumed very quickly and further 
advancement of the dissolution front is only possible when one of the pores increases its flow rate at the expense of 
its neighbors. Such a competition between different flow paths is a characteristic feature of the dissolution problems 
[iMsl, EqI, Ell: a path which gets a slightly larger flow dissolves faster than its neighbors, which decreases its 
hydraulic resistance and makes the flow there even higher. As a result, an intricate, fractal-like structure of channels 
is formed, where all the flow and dissolution is focused. The structure becomes increasingly more branched as Dag// 
is increased, with individual branches again competing with each other for the available flow, which leads to the final 
appearance of a single flowpath joining the inlet and outlet. 

The G parameter affects the dissolution process in a more subtle way, much less evident compared to the effects 
caused by Dag//. For larger values of G we observe more branching in the wormholes and less broadening around 
the inlet, particularly at higher effective Damkohler numbers. This behaviour is connected with a simple mechanism 
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Figure 7: Phase diagram of dissolution patterns for different Dae// and G in the system with hne inlets and outlets. The 
simulations were performed under the fixed flow conditions. The frames in the figure correspond to either the breakthrough 
moment or to the moment when one of the channels has increased its initial diameter 300 times (the frames in a lower left 
corner). The corresponding values of the dimensionless time (|19p are shown for each frame. 

that leads to a drastic increase of dissolution rate in the widest pores in the reaction {G 0) limit - the wider the 
pore is, the larger its surface area {S = iidl) becomes, so that it consumes more reactant compared to other pores. 
This leads to a situation in which increasingly larger amount of reactant is consumed in the first few pores in a given 
flowpath, and the reactant withdraws from the more remote pores (see also a more detailed discussion of this effect 
in the Appendix). This effect is most pronounced when a constant flow is imposed on the system. In the constant 
pressure case, the flow rate increases during the dissolution, hence Dag// decreases and the reactant is penetrating 
deeper inside the system, not being spent in the first few pores. 

This effect is strongly reduced at larger values of G, where the dissolution relatively quickly switches to a transport- 
limited mechanism, which tends to slow down the dissolution rate as the pore opens, due to the inverse dependence 
of keff on the diameter, as observed in Eq. (j9j). 

The geometry of Fig. O with point inputs and outputs, is relevant to a number of karst formation problems 
[26l . [27I . l42j . The point inputs then represents particular sinks through which the water may enter the bedding 
plane, whereas outputs correspond to the springs draining the area. Ewers [26| has also modelled such a situation 
experimentally by creating an artificial model of a bedding plane between the soluble plaster block and a transparent 
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Figure 8: The distribution of wormhole lengths, n{Luj), at the breakthrough time for Dae// = 0.5 and G = 0.1. The soUd Une 
corresponds to the fit n{Lw) ^ . Wormholes shorter than 0.03L and longer than 0.75L (with L the system length) were 
not taken into account in the distribution. The lengths are normalized by the system length and n is normalized by the total 
number of wormholes in the sample, ntot. Results correspond to a simulation with fixed total flow through the system. 



insoluble lower boundary and flushing the system with water through a number of point inlets and outlets. An 
example result of such an experiment is presented in Fig. [6] and bears a lot of resemblance with the patterns obtained 
by use of a network model 

Another experimentally important geometry is the one in which the reagent enters through the entire top surface 
and leaves through the bottom. This is the geometry relevant to the experiments and simulations by Hoefner, Fogler 
and Fredd [qI, [io| on the dissolution of limestone blocks and also Golfier [sf on the salt dissolution in a Hele-Shaw cell. 
In this geometry there are no preferred flow paths along which water may flow more rapidly than through others, 
thus initially many independent wormholes are formed which then compete for the available flow. 

Results of numerical simulation for different values of Dag/ / and G for this geometry are presented in Figure [71 
Constant total flow through the system has been imposed in these simulations (the development of the patterns in the 
analogous simulations under constant pressure conditions is represented in the movies in the Supplementary Material 
[43}). Even though an identical initial geometry has been used in each of the runs, one observes that position of 
the main wormhole varies for different Dag// and G, illustrating sensitivity of the dissolution process to physical 
parameters beyond the geometry. 

The overall dependence of the dissolution patterns on Dag/ / and G shares many similarities with that observed for 
point inlets and outlets. Again, the effective Damkohler number is the key parameter controlling the penetration length 
- for small values of Dag//, as the penetration length becomes comparable with the size of the system, dissolution is 
uniform. Then, as Dag// is increased, a well-defined reaction front appears, which becomes unstable in the course 
of evolution. The linear stability analysis shows that the wavelength of this instability decreases with increasing 
Damkohler number [l^, which is consistent with the results of present simulations - for Dag// = 10~^ one observes 
long- wavelength undulations of the dissolution front whereas for larger Dag/ / separated wormholes are formed where 
the majority of the flow and reaction is focused, while most of the pore space is bypassed. To further investigate 
the scaling of the instability wavelength with the Damkohler number we conducted a series of dissolution simulations 
of regular lattices with uniform pore lengths Iq but with a very small random noise in the initial pore diameter 
distribution: Ad/ do = 10~^. The results show sinusoidal modes developing in the dissolution front, as illustrated 
in Fig. [TOl The wavelength is inversely proportional to the Damkohler number, confirming that penetration length, 
Ip = ^o/Dag//, is the only important length scale in the early stages of dissolution, which is a common feature of the 
convection-dominated dissolution processes [3, 113 • 

The wormholes again compete in a process in which the longer wormholes drain flow from the shorter ones, limiting 
their growth. This time however, in contrast to the point inlet case, due to the translational invariance of the 
system the process of wormhole competition becomes self-similar. The characteristic length between active (growing) 
wormholes increases with time, while the number of active wormholes decreases, which leads to a scale-invariant, 
power-law distribution of wormhole lengths, 

n(L™) ^ L-", (24) 

were n{L^) denotes number of wormholes longer than L^. The fits to Eq. were performed in the range 0.5 > 
Dag// > 5 where the wormholes are well-pronounced and there are relatively many of them, which allows for the 
adequate statistics. Wormholes shorter than 0.03L and longer than 0.75L (with L the system size) were not taken 
into account in the distribution: the former because their lengths are influenced by the lattice discretization effects. 
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Dae//\G 


0.1 


1 


10 


0.5 


1.01 ±0.02 


0.98 ± 0.03 


0.96 ± 0.05 


1 


0.98 ±0.03 


0.90 ± 0.02 


0.89 ± 0.02 


2 


1.01 ±0.02 


0.91 ±0.03 


0.88 ± 0.02 



Table I: Parameter a for different Dae// and G. 




Figure 9: Same as in Fig. [71 but with pore merging (with pore aspect ratio do/lo = 0.1). The frames in the figure correspond 
to the breakthrough moment except for the two frames in the lower left corner, where, for the sake of legibility, we stop the 
simulation at the moment when half of the system has dissolved. 

the latter - because the longest wormholes remain active and the selection process there has not yet been concluded. 

The values of the exponent a obtained from the fitting procedure over a range of different Dag/ / and G numbers are 
presented in Table T The exponent vary slightly with G with the largest values corresponding to the reaction limited 
case (small G). An example distribution of the wormhole lengths obtained at Dag// =0.5 and G = 0.1 is presented 
in Fig. [8] together with the respective fit. The reported values of a in reaction-limited regime (small G) are in a good 
agreement with the 2d Laplacian needle growth model of wormhole- wormhole competition [45|, which yields a = 1. 

Finally, let us discuss the influence of pore merging on the dissolution patterns. The data in Figs. [5] and [71 have been 
obtained in the simulation which has not allowed for the pore merging. For comparison, Fig. [9l shows the patterns 
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obtained using pore merging simulation (with the initial pore aspect ratio d^/l^ = 0.1). Comparing Figs. [7] and [9] 
one notes that merging affects mainly the patterns corresponding to the large values of Dag//. For small Dag// the 
dissolution is uniform, and merging takes place mostly behind the scalloped dissolution front, hence it is not affecting 
the front instability dynamics. Then, as Dag// increases, the merging intensifies, in particular near the inlet and 
along the wormholes. Finally, for large Dag// and small G the non- merging and merging patterns are contrastingly 
different - whereas in the non- merging case dissolution concentrates at the inlet, merging allows the matrix at the 
inlet to become completely dissolved and a steadily advancing front appears, separating fully dissolved pore space 
from the undissolved one. 
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Figure 10: Concentration profiles in a dissolving pore network with a small initial noise, 5d/dQ = 10 ^. The shades of gray 
represent the magnitude of the concentration field with black standing for c = Cin and white for c < O.OTcin- 



IV. DISSOLUTION OF THE REGULAR LATTICE 

In this section we consider the dissolution of an ideal hexagonal network. Beside being an interesting problem of its 
own, it will also give us insight into the mechanisms governing the formation of the dissolution patterns. Let us then 
set initially all the diameters and wormhole lengths in the network to uniform values, d = do and / = Iq throughout 
the whole system. Then, to induce a localized growth of such a system, we make a single small cut in the inlet region, 
increasing the diameter by a factor of four in a vertical column comprising 10 network nodes situated at the center of 
the inlet manifold (cf. Fig. [11]). Fig. [11] presents the dissolution patterns obtained for G = 1 and different Dag// in 
such a system (for the simulation without pore merging). Two strikingly different patterns can be observed there: for 
small Dag// a dendrite- like wormhole is formed, with a characteristic, spearhead-like shape. In this case the width 
of the wormhole is much larger than the pore scale. Around Da^// = 0.54, a sudden transition takes place, to an 
'inverted Y' structure, involving two channels growing sideways at equal angles from the main branch. The channels 
are only one pore wide and almost all the dissolution and flow becomes quickly concentrated there. The value of G 
affects the patterns to a weak extent only, its main effect being to shift the transition point between the patterns: for 
G = 10 the critical Da^ff value becomes Da^// = 0.67, whereas for G = 0.1 it is Da^// = 0.35. This is in agreement 
with the observations discussed in Sec. HIT] larger G values Hmit the dissolution at the inlet and thus increase the 
react ant penetration length and allow for the formation of more diffuse dissolution patterns. 

To elucidate the origin of the patterns observed in the dissolution of regular lattices, let us consider a simplified 
triangular network shown in Fig. [12] Here A represents the tip of the cut, with pA = ^ whereas D represents the outlet 
boundary of the system, with pn = 0. By solving the Kirchoff equations in this case, one gets p|Zp^ = |- The flow 
along side pores is therefore slightly larger than in the central ones at the beginning of the simulation, when resistances 
of all the pores are the same. However, at large Dag//, when almost all the reactant is spent in the first pores, even a 
small difference in the flows will be enhanced by the dissolution. As a result, the hydrodynamic resistance of a channel 
BC will decrease much faster than that of BX, and the concentration of reactant at G will be significantly larger than 
that in X. Thus the line AD will continue to grow and soon all of the flow in the system will be concentrated there. 
On the other hand, for smaller Dag// both BX and BC will dissolve almost uniformly. In that case, the concentration 
of reactant in X will be in fact larger than that in G (since X is being fed by two reactant-bearing pores instead 
of one). The dissolution will then be concentrated in the central channels, and a dendrite-shaped wormhole will be 
formed. 

We expect the shape of the dendrite in this case to be largely independent of the underlying lattice, and in fact 
similar shapes are observed in the case of random lattices in Fig. [7] For a regular lattice, an analytical derivation 
of such a shape should in principle be possible, which remains a future task. On the other hand, the inverted Y 
pattern formed at large Dag// is strongly connected with the underlying network. The arms of the 'inverted Y' figure 
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Figure 11: Dissolution patterns of a regular network for different Dae// and G = 1. Initially, a short channel was cut in the 
center of inlet manifold to induce the localization of the flow. 




Figure 12: A simplified triangular network with a single input node (A) and line output (D). 

propagate along the lattice directions and the pattern will look differently if the network is changed. It is also strongly 
unstable, as the arms begin to compete with each other, as can be observed in Fig. [TTJ 

The above considerations, although originally pertaining to the regular network, have, in fact, a more general 
relevance, since they illustrate two basic mechanisms that govern the emergence of dissolution patterns in the network. 
At small or intermediate Dag// diffuse, multi-pore patterns are formed in the extended region where the initial 
concentration of reactant is large. For large effective Damkohler numbers the local pressure drop effects become more 
important and the dissolution is focused along a thin pore-wide wormhole which is more irregular and branched, 
reflecting the disorder of initial network. 

V. OPTIMAL INJECTION RATE 

Understanding of the emergence of the dissolution patterns is particularly important in optimization of carbonate 
reservoir stimulation treatments, where the relevant question is how to get the maximum increase of permeability for a 
given amount of reactive fluid. Numerical and experimental investigations of reactive flows in porous media [sl, [lo|, [lO, 
[2l|,[28| suggest that there exists an optimum injection rate, which maximizes the permeability gain for a given amount 
of reactant. If the injection rate is relatively small, a large portion of the reactant is wasted by the unproductive 
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Figure 13: (Color online) Pore volume to breakthrough for different Dae// and G for the simulation with fixed total flow on 
200 X 200 lattice (broadening factor /3 = 4). 

dissolution of the inlet pores while the overall increase in permeability remains small. On the other hand, for very 
large injection rates, the reactant is exhausted on a uniform opening of all of the pores in the system, which is also 
inefficient in terms of permeability increase. The optimum flow rate must give rise to spontaneous channeling, since 
the reactant is then used to create a small number of highly permeable wormholes, which transport the flow most 
efficiently. To quantify the optimization with respect to Dag// and G, we measured the total volume of reactive fluid, 
Vb^ that must be injected into the pores to obtain the breakthrough, which we identify with the moment when the 
diameter of at least one of the outlet pores increases by a factor of /3. It is convenient to express Vb in terms of a 
dimensionless quantity 

where 7 is the acid capacity number ([20]) . Vb is the initial total fluid volume in the network, Q is a (constant) injection 
rate and T5 is the breakthrough time. A diagram of Vb as a function of Dag// and G is presented in Figure [T3l In 
these simulations we did not allow for the pore merging. 

As elucidated above, for both small and large effective Damkohler numbers Vb is relatively large. However, the value 
of G also has a non-trivial impact on Vb. Namely, for small values of G (in reaction-limited regime) the dissolution in 
the inlet region is strongly enhanced, particularly for large Damkohler numbers, which leads to the high consumption 
of reactant there and the associated increase in V5. 

In the real physical system, the value of G cannot be varied during the experiment, since both the diffusion constant 
and reaction rate are material properties. Thus, changing the injection rate moves the system along a line of constant 
G, sweeping the range of Dag// values. This leads to V5(Dae//) dependence, the example of which, for G = 1, is 
presented in Fig. [HI Two different conditions for breakthrough are compared here, corresponding to /3 = 4 and 
/3 = 8, again without pore merging. A stronger breakthrough condition results in an increase of V5, especially for 
smaller Dag//. For larger Dag//, however, a pronounced wormhole is formed in the system and at the moment of its 
breakthrough the dissolution in the outlet pores becomes so dynamic that Vb becomes almost independent of /3. 

The most striking feature of Fig. [14] is the presence of multiple minima on V5(Dag//) curve. To interpret this 
finding, let us first analyse the dissolution of a single, long channel of length L spanning the whole length of the 
system and comprised of = L/Iq elementary pores. As derived in the Appendix, if initially the channel is uniform, 
d{x,t = 0) = do and the concentration at the inlet is given by c{x = 0,t) = Qn, then the diameter of the channel 
d{x, t) is given implicitly by (|A13p 

(1 + G)Da^^r^- + Gin I | - 2 farccoth(d) - ^icco\h{din(t))] = 0, (25) 

where d = d/do is the diameter of the channel scaled by the initial diameter do and din{i) is the diameter at the inlet 
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Figure 14: Breakthrough volume as a function of Dae// for G — 1 for the simulation with fixed total flow on 200 x 200 lattice. 
The data for two different conditions for breakthrough is presented: /3 = 8 (circles) and /3 = 4 (squares). The lines are a guide 
to the eye. 



{x = 0) which evolves according to (jAlSp 

a/2G£+G(G + 2) + 1-1 

d^n = ^ ^ . (26) 

Finally, the effective Damkohler number for the channel is given by the formula analogous to (|2T]), i.e. 

Defining breakthrough time T5 as before as the moment at which the end of a channel {x = L) has broadened /3 
times, we obtain V5(Da^^}^^) curves presented in Figure [15] for different G regimes (G = 0.1, 1 and 10). Interestingly, 
even on a level of a single tube we observe the minimum in H(Da:;^^"") dependence. The increase of Vb at large 

Da^jj^^ is caused by the nonuniform dissolution of the tube in the low- flow limit: the inlet dissolves faster than the 
outlet, which leads to the increase of the reactive surface area in the inlet region and further depletion of the react ant 
from the downstream regions. On the other hand, the increase of T4(Da^^;^"") at very high flows (low Da^jj^^) is 
connected with the fact that a significant portion of the reactant is then simply flushed through the system, not 
reacting with the walls. The minimum in V^iDdiffJ^'^) is present for all values of G, both in reactive-limited and 
transport-hmited regime, though its position shifts to larger Da^jj^^ as G is increasing. 

It is interesting to note that although in general (|25]) cannot be solved explicitly, for a number of specific values 
of G a closed form expression for T5(Dae//) can be found. These include, among others, G = and G = 00, i.e. 
reaction-limited and transport-limited case respectively, which are worked out in the Appendix, but also G = 1, which 
corresponds to the case when both diffusion and reaction are important. Putting x = L and d = j3 into ([25]) we obtain 
in the latter case 

T5 = ^y^e^^^/r (4 + (/3 - l)e''<yr) (28) 
or, in terms of the volume of reactive fluid injected, V5 = qchannTb = qchanndoTb/2kj 

Vb = ^= e^-e?r(4 + (/3 - l)e^^ejr^), (29) 

Vo 2Da^)p y J y J 

where we have normalized V5 by the initial volume of the channel, Vb = irLdl/^. The above corresponds to the solid 
line in Fig. [15] 

The single-channel model should adequately describe the results of the network dissolution in the regime of large 
Dag//, where the competition between the flowpaths is strong and leads to the emergence of a solitary winning 
wormhole. Figure [16] shows the profile of a channel formed in such a regime (for Dag// = 5 and G = 1) at the 
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Figure 15: Pore volume to breakthrough as a function of Da^^^^^ for C = 0.1 (reaction-hmited case, dashed hne) G — 1 
(mixed case, sohd) and G = 10 (diffusion-hmited case, dot-dashed) based on the analytical solution for the profile of a single 
dissolving channel (|25)) . The value of the broadening factor used here is /3 = 4. 




Figure 16: Dissolution of a channel obtained by numerical simulation on the network (empty circles) and exact analytical 
solution for single wormhole, Eq. ([25]) (dashed line). Comparison was made for Dae// = 5 and C = 1 at the breakthrough time. 



breakthrough time. As observed, the profile of the channel agrees with the analytical result (|25]) very well, except in 
the region of large x. This is caused by the fact that in the vicinity of channel tip there is an intensive leakoff from the 
channel towards the matrix [40|, thus the assumption of the constant flow throughout the channel, used in derivation 
of ([25]) ceases to be valid. 

At smaller Damkohler numbers, there are many alternative flowpaths between inlets and outlets which divide the 
available flow between themselves. In the simplest approximation, we can assume that the flow is divided equally 
between the flowpaths (c/. Fig. [TT]) and then again apply Eq. ([25]) to one of such flowpaths. However, in order to 
calculate the Damkohler number ([27]) in one of these channels we need to assess what percentage of the total flow 
goes through a particular path. This depends on the amount of focusing in the system, which, as demonstrated in 
Sec. IIIH is a function of the overall Damkohler number. Dag//, characterizing the whole system: for small Dag// the 
dissolution is uniform throughout the system and Qchann will be equal to Q/Niniet- For larger Dag//, undulations 
begin to be formed on the dissolution front, but the emerging fingers are relatively thick, comprised of many parallel 
pores. For even larger Dag// the fingers become narrower and finally, for the largest Dag//, only one pore at each 
height carries all the fluid, and qchann = Q- 

In general, the percentage of the flow focused in the active flowpaths, qchann/ is an increasing function of Dag// 
of a rather complicated nature. Fig. [181 presents an example of such a function numerically estimated for G = 1 
together with a phenomenological flt to the data: 

/I (lnDa^^^+a3)2 

{qchann /Q) = \n {N inlet) (^-(tanh(ai + In Dag//) - 1) - -4 j, (30) 

with ai 4.10, a2 0.42, as 2.34, and 3.03. 

An interesting feature of the dependence <7c/iann (Dag//) is a plateau for intermediate values of Dag// (starting from 
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Figure 17: A schematic of a simple model of disconnected channels. 
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Figure 18: The characteristic flow through an active flowpath, Qchann, as a function of Dae// together with a phenomenological 
fit, Eq. (|3Q)). 



Dag// ^ 0.02). The analysis of the patterns in this region shows that it corresponds to the onset of branching of the 
dissolution channels. However, both the number of active channels and their widths remain approximately constant 
up to about Dag// ~ 0.2. 

Combining the Qchanni^^eff) dependence of Eq. (|3Q|) with (|28|) leads to 



3-1 ^-effN.Q JD^effN.Q 



(31) 



where we have linked the channel Damkohler number (|27|) to the network Damkohler number ([16]) by 



chann 
eff 



nkdnL 



Qchann (l + G) qoil + G) Qcha 



2Nyq 

chann 



(32) 



In the above, we have used the fact that L = IqNx and Q = NinietQo = ^TV^go^ since there are two channels leading 
out of each inlet node (c/. Fig. [2]). In terms of the injected volume (still for G = 1) 



4i^har^r^ (4 + (/3 - l) e^''y^ohar^r^ ) ^ 



(33) 



where V5 is normalized by the initial volume of the pores in a hexagonal x Ny lattice, Vq = '^NxNyirlod^/ A. Note 
that Qchann the above formulae is in itself a function of both Da^ff and the system size, as given by Eq. (|3Q|) . 

The resulting Vi){Daeff) dependence for 200 x 200 system is presented in Fig. [191 Additionally, we plot here 
the V5(Dae//) dependences obtained under the assumption that the entire flow in the system is always focused in 
a single channel, irrespectively of the Damkohler number, i.e. Qchann = Q and yet another curve corresponding to 
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Qchann = Q/^iniet^ i-G. with the assumptioH that there is no flow focusing at ah and the fluid flows through ah 
available flowpaths. The comparison of these dependences helps us to understand the origin of particular minima in 
V5(Dae//) curve. The deepest minimum at Dag// ^ 1 corresponds to the optimal condition of dissolution in a single 
channel model with Qchann = i-e. to the case when the entire flow is focused in a single flowpath. On the other 
hand, the shoulder at the left, near Dae// ^ 0.002 corresponds to the optimal conditions of dissolution in the case 
when all the flowpaths dissolve uniformly. Finally, the local minimum at Dag// ^ 0.03 corresponds to the plateau 
region on ^c^ann (Dag//) curve and is an effect of superimposing the dependences in Figs. [T5l and [TSl All of the 
above-discussed features are also visible in the V5(Dae//) curves constructed from the numerical data in FiglMl 




10-4 0.001 0.01 0.1 1 10 

D^eff 



Figure 19: Pore volume to breakthrough obtained for a simple model of non-interacting channels, using Eq. (|25|) with the 
flow in the active channel, Qchann obtained based on Eq. (|30p . The dotted line corresponds to (|25p with Qchann = Q (all flow 
focused in a single channel) whereas the dashed line corresponds to the uniform dissolution with Qchann — Ql^iniet 



The data discussed so far has been obtained without pore merging. Since the larger channel grows faster than 
two smaller ones, merging will in general decrease the breakthrough times. Fig. [20] presents the breakthrough curves 
for the simulation with pore merging for two different values of initial aspect ratio, d^jl^ = 0.1 and d^jl^ = 0.025. 
Comparing these curves with that for a non- merging case, one observes a signiflcant speed-up of dissolution and 
decrease of breakthrough times in the intermediate range of Dag//. For small Damkohler numbers the impact of 
merging on the dynamics is much weaker, since the pores then remain thin and thus hardly merge. On the other 
hand, for large Dag//, there emerges a single wormhole which eats up the neighboring ones. In the merging simulations 
the main wormhole is broadening faster than without merging, particularly near the inlet. As already mentioned, this 
leads to a large consumption of reactant there which slows down the advancement of the dissolution front and gives 
rise to the increase of breakthrough times. 

The main difference between the breakthrough curves in merging and non-merging case is however a gradual 
disappearance of multiple minima structure. Already at (io/^o = 0.025 the barrier between the two main minima 
becomes lower and at do/^o = 0.1 it vanishes altogether and a single minimum at Dag// — 0.2 is observed. A similar 
dependence of V5 on Damkohler number has been reported previously [sl, [lo|, [20|, HH, [28| , both in the experiments and 
simulations of porous media dissolution. 

Finally, we discuss the effects of the system size on the above results. Fig. [211 shows the Vd{Daeff) curves for the 
simulation without pore merging for various N x N lattices with TV ranging from 20 to 400. We flnd it convenient to 
plot NVd as the dependent variable here, as then the curves asymptote to a single curve as Daeff 0. As observed, 
for small system sizes only one minimum appears in Vd{Daeff) dependence. For larger sizes this minimum bifurcates 
into two minima, and then, as the size is further increased, a third shallow minimum appears for small N Daeff values. 
This behavior can be further elucidated using the semi-analytical model of Eqs. (j25ti33p . as presented in Fig. [22l In 
particular, the Dag// limit of ([3T]) gives a small Daeff asymptote in the form 

(/3-l)(3 + /3) , , 

NV, ^ ^-^^—^ ^, (34) 

3i;ag// 

indeed independent of N. Moreover, as observed in Fig. [22] in the large Daeff limit the curves also converge, this 
time to 

NVb ^ y ~^^' e^^-^^ (35) 



3Da, 



eff 
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Figure 20: Breakthrough volume as a function of Dae// for G — 1 and f5 — A for the simulation without pore merging (circles) 
and with pore merging. In the latter case, two values of the initial pore aspect ratio are considered: do/lo — 0.025 (triangles) 
and (Iq/Iq — 0.1 (squares). The simulations were carried out on 200 x 200 lattice. 

In the numerical simulations of Fig. [21] this limit is hard to attain, since at the large Damkohler numbers the dissolution 
is highly nonuniform which is a source of convergence problems in the linear solver. The analytical model shows a 
similar minimum structure to the numerical data - at small N a single minimum is present, at Daeff ~ 1, then the 
second one appears, at Daeff ^ 0.03 and finally the third at Daeff ^ 0.001. As before, the rightmost minimum 
corresponds to the case when the entire flow is focused in a single wormhole, the central minimum corresponds to the 
plateau region on ^c/iann (Dag//) and the leftmost minimum is connected with a uniform dissolution of the network. 
Importantly, the positions of these minima are almost entirely independent of the size of the system. 

It is worth to note, however, that the model of Eqs. (|30ti33|) gives only a qualitative agreement with the numerical 
data - even though the general three-minima structure is reproduced correctly, there are significant differences in 
relative depths of the minima and the heights of the barriers separating them. The main approximation in the model 
is the assumption that the percentage of the flow focused in the main channel (|3Q|) remains constant throughout the 
simulation. In reality, however, initially all the channels have similar flow rates, and then gradually the flow focuses 
in the few main channels. Thus the model of Eqs. (jSOtiSSp will in general overestimate the total dissolution rate 
in the main channels throughout the simulations, particularly for large Damkohler numbers, where flow focusing is 
the strongest. This leads to much larger depth of the rightmost minimum in T4(^<^e//) than that observed in the 
simulations. 

Finally, let us discuss the size effects for the pore-merging simulations. The data presented in Fig. [23] shows that 
in this case there are only slight differences in Vh{Daeff) dependences as N is increased. All of the curves have a 
single minimum only, which seems to be drifting towards smaller Daeff values as the size is increasing, however the 
minimum is rather wide so that in the entire region 0.1 < Daeff ^ 1 the values of breakthrough volumes are similar. 

VI. SUMMARY 

In this paper we have studied dissolution of porous medium using an evolving network model. The relevant 
dimensionless parameters to characterize dissolution patterns in this system are the effective Damkohler number. 
Dag//, relating the reaction rate to the mean fluid velocity in the pores and G, which measures the relative importance 
of diffusive and reactive phenomena. As Dag// is increased, a transition from uniform dissolution to the strongly 
localized flow is observed, with the wormholes appearing at about Dag// ^ 0.02. For intermediate values of Dag// 
(0.02-0.5) the wormholes are diffuse and involve many parallel pores. However, as Dag// is further increased, the 
wormhole diameter is drastically decreased and dissolution becomes focused along thin, pore-wide, highly branched 
fiowpaths. This transition in wormhole shape is particularly dramatic and abrupt in the case of regular network, 
where a characteristic, 'inverse Y-shaped' channels appear at Dag// ^ 0.5 (for G = 1). 

We have also analyzed the problem of finding an optimum flow rate that gives a maximum increase in permeability 
for a given amount of react ant. The dependence of Vb (the volume of the fluid needed to open the system) on Dag// 
turned out to be nontrivial, with a number of local minima, corresponding to different regimes of the dissolution. 
Interestingly, if the pores are allowed to merge together as they dissolve, the Vb{Daeff) curve becomes simpler, with 
a single minimum. 
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Figure 21: Breakthrough volume as a function of Dae// for G = 1 and /3 = 4 for the simulation without pore merging for 
N X N lattice with N = 20 (solid), 100 (dashed), 200 (dotted) and 400 (dot-dashed). 
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Figure 22: Pore volume to breakthrough obtained using a model of non-interacting channels (Eqs. ([25|) - (|3Q|) for x lattice 
with iV = 20 (solid), 100 (dashed), 200 (dotted) and 400 (dot-dashed). 




Figure 23: Breakthrough volume as a function of Dae// for G = 1 and /3 = 4 for the simulation with pore merging {do/lo = 0.1) 
ioT N X N lattice with N = 20 (solid), 100 (dashed), 200 (dotted) and 400 (dot-dashed). 
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Appendix A: The growth of a single cylindrical channel 

Let us consider the problem of dissolution of a single, long channel of length L. The shape of the channel is described 
by its diameter d{x^ t) as a function of axial coordinate x and time t. In principle, one can imagine such a channel as 
a collection of A^^^ = L/Iq pores from Sec. HIl connected in a serial manner. Let us assume that initially the channel is 
uniform, d{x^ t = 0) = do and the concentration at the inlet is given by c{x = 0, t) = q^- 

The erosion equation, (p!3|) can then be written as 

Otd=- -c (Al) 

where 7 = Cin/h'CsoU whereas G = and the dimensionless fields c = c/cin^ d = d/do has been introduced. Next, 
the concentration balance. Eq. (pTj) gives 

q l^Gd ^ ^ 



Scahng time by dQ/2kj and the axial coordinate by the channel length: 

2k2 
do 



(A3) 



L 

leads to 



x=^ (A4) 



and 



dfd = — —.c (A5) 
' l^Gd 



g^c = = - (1 + G)Da^^p^^,c, (A6) 



q l^Gd l^Gd 
where the effective Damkohler number for the channel is given by 



j>.^chann _ T^dokL / Qchann / * 

^^eff - YTG ^ ^ 



Differentiating ()A5p with respect to x and using ()A6p we get 



%d = -(1 + G)Da^f."" ^-^c ^—^cdj (A8) 



Finally, expressing c in terms of d^d using again (jA5j) leads to the equation in terms of d only 



d^fd = -(1 + G)Dag^r^ . dtd —.dfddj (A9) 

-^-^ 1 + l^Gd 



or, regrouping the terms 



d^t \^d{l + G^)j + dt Q(l + G)Da^^7"d2 j ^ (AlO) 
Integrating over time leads to the relation 

d, (d{l + G^U + i(l + G)Da^5p(P ^ (All) 
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The function W{x) can be obtained from the initial condition: d{i = 0) = 1, which gives 



2 

The solution of (jAlip can only be obtained in an implicit form: 

' d2 _ 1 ^ 



W = 1(1 + G)I)al'}'}"" (A12) 



(1 + G)Da^^;;""i + G In ( J - 2 (arccoth(d(x, i)) - arccoth(di„(t))) = (A13) 



where din{i) is the inlet diameter, which is the solution of ()A5p with c = 1, i.e. 



Sid™ = ^ (A14) 

1 + Gdin 



which gives 



/2Gt + G(G + 2) + l-l 
d^i) = ^ (A15) 

We shall now consider two limiting cases of Eq. (|A13[) : the reaction-limited case, G ^ 0, and the transport Hmited 
case, In the reaction limit ()A13p takes the form 

j^^chann^ _ ^ |^arccoth(d) - arccoth(din)) = (A16) 

where 

j^^chann ^ ^-^ Da^^^^ = ^^^^ (A17) 

G^o q 

In this case ()A15p becomes simply 

din{i) = l^i (A18) 
The above two equations can then be combined to yield 

d = coth Qoa^^^^^x + arccoth(l + t)^ (A19) 

Interestingly, the above profile has a well-defined limit as t ^ oo, namely 

^^^oo = coth Qoa^^^^^x^ . (A20) 

Thus, at each given point, the channel ceases to grow after a certain time. This reflects the fact that as the diameter of 
the channel increases, there is a growing reactive surface, and hence the concentration is consumed faster, particularly 
near the inlet. Thus, paradoxically, the concentration profile should recede as time goes on. The receding concentration 
profile can indeed be observed, since in the reaction-limited case 

c = dtd = (^(i + 1) sinh Qoa^^^^^f ^ + cosh Qoa^^^^^f ^ ^ (A21) 

which is a decreasing function of time. 

In this limit, we can also calculate analytically the breakthrough time, i.e. the time such that d{x = 1, f = T5) = {3. 
Inverting Eq. ()A19p leads to 

n = coth rarccoth(/3) - ^Da"^^^^^ - 1 (A22) 
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which can be rewritten as 

fl + cothf|Da"^^^^)) 

n = 4^ -^^ ^ (A23) 

coth UDa"^^^^ j - (3 

The above formula can be apphed only if the denominator is positive, which gives the condition 

j^^chann ^ 2arccoth(/3) (A24) 

a direct consequence of an upper limit of diameter growth, as given by (|A2Q|) . 
In the transport limited case, G ^ oo, Eq. ()A13p takes the form 



d=\jl^ gg-TTDShLf/g (A25) 

where we have used the fact that limc^oo Da^jj^^ = TiD?>\iL/q. This time the diameters increase without a limit. 
The breakthrough time is then given by 

= 1g {(3^ - 1) e-^ShL/g ^^26) 
This time T5 is well-defined in the entire fiow range, but it dramatically increases as g — 0. 
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